cohortfile <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/cohortfiles/dti_fa/dti_fa_cohortfile.csv")
cohortfile <- cohortfile %>% rename(subjectID=sub) %>% select(subjectID, age, sex, mean_fd)
tract_profiles <- fread("/cbica/projects/luo_wm_dev/input/HCPD/HCPD_tractprofiles/all_subjects/collated_tract_profiles_reoriented.tsv")
cohortfile <- cohortfile %>% select(subjectID, age, sex, mean_fd)
gam_df <- merge(tract_profiles, cohortfile, by="subjectID")
gam_df <- gam_df %>% mutate(tract_node = paste0(tract_hemi, "_", nodeID))
gam_df$sex <- as.factor(gam_df$sex)subjects <- read.table("/cbica/projects/luo_wm_dev/input/HCPD/subject_list/HCPD_subject_list_N569.txt")
subjects <- subjects %>% setNames("subjectID")
subject_subset <- subjects$subjectID[c(1:75)]
gam_df <- gam_df %>% filter(subjectID %in% subject_subset) # fit GAMs on a subset of 75 participantsfit_tractwise_gamm_re <- function(tract_name, scalar) {
df <- gam_df %>% filter(tract_hemi == tract_name)
df$subjectID <- as.factor(df$subjectID)
te_gam <- gamm(as.formula(sprintf("%s ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') +
ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd", scalar)), random = list(subjectID=~1), data = df)
print(summary(te_gam$gam))
print(k.check(te_gam$gam))
print(paste(tract_name, "done"))
return(te_gam) # return the model object
}
fit_tractwise_gam_re <- function(tract_name, scalar) {
df <- gam_df %>% filter(tract_hemi == tract_name)
df$subjectID <- as.factor(df$subjectID)
te_gam <- gam(as.formula(sprintf("%s ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') +
ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd + s(subjectID, bs = 're')", scalar)), data = df)
print(summary(te_gam))
print(k.check(te_gam))
print(paste(tract_name, "done"))
return(te_gam) # return the model object
}
fit_tractwise_bam <- function(tract_name, scalar) {
df <- gam_df %>% filter(tract_hemi == tract_name)
df$subjectID <- as.factor(df$subjectID)
te_gam <- bam(as.formula(sprintf("%s ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') +
ti(age, nodeID, k=c(3, 24), fx = F) + s(nodeID, by = subjectID) + sex + mean_fd", scalar)), data = df, method = "fREML", discrete = TRUE)
print(summary(te_gam))
print(k.check(te_gam))
print(paste(tract_name, "done"))
return(te_gam) # return the model object
}tractwise_gamm_re <- lapply(unique(gam_df$tract_hemi), fit_tractwise_gamm_re, scalar="dti_md")
names(tractwise_gamm_re) <- unique(gam_df$tract_hemi)
saveRDS(tractwise_gamm_re, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gamm_re.RData")
tractwise_gam_re <- lapply(unique(gam_df$tract_hemi), fit_tractwise_gam_re, scalar="dti_md")
names(tractwise_gam_re) <- unique(gam_df$tract_hemi)
saveRDS(tractwise_gam_re, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gam_re.RData")
tractwise_bam <- lapply(unique(gam_df$tract_hemi), fit_tractwise_bam, scalar="dti_md")
names(tractwise_bam) <- unique(gam_df$tract_hemi)
saveRDS(tractwise_bam, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_bam.RData")tractwise_gamm_re <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gamm_re.RData")
tractwise_gam_re <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gam_re.RData")
tractwise_bam <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gam_smooth.RData")I fit 3 tractwise models on a subset of HCP-D (full sample = 569; subsetted sample = 75 that represents the full age range)
Using gam with subjectID as a random effect:
gam(dti_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd + s(subjectID, bs = 're'), data = df)
Using gamm and having subjectID as a random effect (gamm relies
on lme, I believe):
gamm(dti_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd, random = list(subjectID=~1), data = df)
Using bam and using separate smooth function for subjectID
bam(dti_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd + s(nodeID, by = subjectID), data = df, method = "fREML", discrete = TRUE)
plot_compare_residuals <- function(tract, num_subjects) {
gam_df_toplot <- gam_df %>% filter(tract_hemi == tract)
# Extract residuals for each model
gam_df_toplot$residuals_gam <- resid(tractwise_gam_re[[tract]])
gam_df_toplot$residuals_gamm <- resid(tractwise_gamm_re[[tract]]$gam)
gam_df_toplot$residuals_bam <- resid(tractwise_bam[[tract]])
subjects_to_plot <- subject_subset[c(1:num_subjects)]
# Plot residuals for each model
gam_plot <- ggplot(gam_df_toplot[gam_df_toplot$subjectID %in% subjects_to_plot, ], aes(x = nodeID, y = residuals_gam, color = factor(subjectID))) +
geom_line() + theme(plot.title = element_text(hjust=0.5, size = 16, face="bold"),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
labs(title = "Raw Residual Curves from Model 1 \ngam: s(subjectID, bs = 're')", x = "Node ID", y = "Residuals")
gamm_plot <- ggplot(gam_df_toplot[gam_df_toplot$subjectID %in% subjects_to_plot, ], aes(x = nodeID, y = residuals_gamm, color = factor(subjectID))) +
geom_line() + theme(plot.title = element_text(hjust=0.5, size = 16, face="bold"),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
labs(title = "Raw Residual Curves from Model 2 \ngamm: random = list(subjectID=~1)", x = "Node ID", y = "Residuals")
# Plot smooth residuals from the new approach
bam_plot <- ggplot(gam_df_toplot[gam_df_toplot$subjectID %in% subjects_to_plot, ], aes(x = nodeID, y = residuals_bam, color = factor(subjectID))) +
geom_line() + theme(plot.title = element_text(hjust=0.5, size = 16, face = "bold"),
axis.title.x = element_blank(),
axis.title.y = element_blank()) +
labs(title = "Smooth Residual Fit from Model 3 \nbam: s(nodeID, by = subjectID)", x = "Node ID", y = "Residuals")
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=16))
y.grob <- textGrob("Residuals",
gp=gpar(col="black", fontsize=16), rot=90)
z.grob <- textGrob(tract,
gp=gpar(col="black", fontsize=20))
plots <- ggarrange(gam_plot, gamm_plot, bam_plot, ncol=3, common.legend=T, legend = "right")
return(grid.arrange(arrangeGrob(plots, left = y.grob, bottom = x.grob, top = z.grob)))
}tracts <- unique(gam_df$tract_hemi)
tracts <- tracts[-c(which(tracts=="Forceps_Major"))]
lapply(tracts, plot_compare_residuals, 5)## [[1]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[2]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[3]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[4]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[5]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[6]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[7]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[8]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[9]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[10]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[11]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[12]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[13]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[14]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[15]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[16]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[17]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[18]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[19]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[20]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[21]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## [[1]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[2]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[3]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[4]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[5]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[6]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[7]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[8]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[9]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[10]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[11]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[12]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[13]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[14]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[15]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[16]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[17]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[18]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[19]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[20]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
##
## [[21]]
## TableGrob (1 x 1) "arrange": 1 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
plot_compare_acf <- function(tract) {
# Define the layout: 2 rows (1 for title, 1 for plots), 3 columns for plots
layout(matrix(c(1, 1, 1, 2, 3, 4), ncol = 3, byrow = TRUE), heights = c(1, 4))
# Plot the title
par(mar = c(0, 0, 1, 0)) # Set margins for the title
plot.new()
text(0.5, 0.5, tract, cex = 2, font = 2)
# Reset margins and set up for 1 row, 3 columns for the plots
par(mar = c(5, 4, 4, 2) + 0.1)
# Plot the ACFs
acf(resid(tractwise_gam_re[[tract]], type = "response"), main = paste0("Model 1 \ngam: s(subjectID, bs = 're')"))
acf(resid(tractwise_gamm_re[[tract]]$gam, type = "response"), main = paste0("Model 2 \ngamm: random = list(subjectID=~1)"))
acf(resid(tractwise_bam[[tract]], type = "response"), main = paste0("Model 3 \nbam: s(nodeID, by = subjectID)"))
}## [[1]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.972 0.934 0.892 0.847 0.802 0.761 0.725 0.696 0.674 0.658 0.647 0.641
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.637 0.635 0.634 0.634 0.633 0.633 0.631 0.629 0.624 0.618 0.610 0.600 0.590
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.579 0.569 0.558 0.549 0.541 0.533 0.526 0.519 0.512 0.504 0.496 0.487 0.478
##
## [[2]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.971 0.936 0.899 0.862 0.825 0.792 0.764 0.741 0.723 0.710 0.701 0.695
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.691 0.687 0.685 0.682 0.679 0.675 0.669 0.662 0.651 0.638 0.624 0.609 0.594
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.579 0.566 0.555 0.545 0.535 0.525 0.515 0.506 0.496 0.486 0.477 0.468 0.459
##
## [[3]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.981 0.961 0.946 0.929 0.910 0.891 0.873 0.857 0.843 0.828 0.815 0.801
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.788 0.776 0.764 0.752 0.740 0.728 0.716 0.705 0.693 0.682 0.671 0.660 0.649
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.638 0.628 0.617 0.606 0.595 0.585 0.573 0.562 0.551 0.539 0.527 0.515 0.504
##
## [[4]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.979 0.958 0.941 0.922 0.901 0.881 0.863 0.846 0.831 0.817 0.803 0.790
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.778 0.765 0.752 0.739 0.726 0.714 0.701 0.689 0.677 0.665 0.653 0.642 0.630
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.618 0.607 0.595 0.584 0.572 0.561 0.549 0.537 0.526 0.515 0.504 0.492 0.481
##
## [[5]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.723 0.719 0.610 0.530 0.494 0.488 0.491 0.485 0.470 0.451 0.432 0.421
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.422 0.433 0.447 0.454 0.452 0.443 0.430 0.418 0.413 0.412 0.415 0.418 0.416
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.410 0.402 0.393 0.385 0.378 0.370 0.360 0.349 0.339 0.330 0.324 0.321 0.320
##
## [[6]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.712 0.693 0.573 0.485 0.450 0.450 0.458 0.455 0.440 0.418 0.401 0.400
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.414 0.438 0.459 0.468 0.463 0.447 0.427 0.414 0.411 0.415 0.421 0.424 0.420
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.411 0.399 0.388 0.378 0.370 0.360 0.346 0.331 0.319 0.312 0.310 0.311 0.313
##
## [[7]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.968 0.917 0.861 0.807 0.761 0.727 0.702 0.685 0.673 0.664 0.659 0.655
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.652 0.649 0.645 0.640 0.632 0.623 0.613 0.603 0.594 0.586 0.578 0.571 0.566
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.560 0.553 0.544 0.533 0.522 0.511 0.502 0.492 0.483 0.474 0.464 0.454 0.445
##
## [[8]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.968 0.916 0.861 0.809 0.764 0.729 0.704 0.685 0.670 0.660 0.653 0.648
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.645 0.642 0.637 0.632 0.625 0.617 0.608 0.599 0.590 0.581 0.573 0.564 0.556
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.548 0.540 0.531 0.523 0.515 0.507 0.500 0.493 0.486 0.478 0.469 0.459 0.448
##
## [[9]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.975 0.944 0.917 0.892 0.865 0.837 0.810 0.784 0.761 0.740 0.722 0.707
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.694 0.682 0.673 0.664 0.656 0.648 0.640 0.634 0.628 0.622 0.616 0.611 0.605
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.600 0.594 0.589 0.583 0.577 0.571 0.563 0.555 0.546 0.536 0.526 0.515 0.504
##
## [[10]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.973 0.940 0.908 0.878 0.847 0.817 0.790 0.767 0.748 0.733 0.721 0.712
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.705 0.697 0.690 0.683 0.675 0.667 0.658 0.648 0.638 0.628 0.617 0.606 0.596
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.587 0.579 0.571 0.563 0.555 0.547 0.539 0.530 0.521 0.512 0.502 0.491 0.480
##
## [[11]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.972 0.946 0.932 0.918 0.900 0.878 0.857 0.837 0.821 0.806 0.793 0.782
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.771 0.761 0.752 0.742 0.733 0.724 0.715 0.706 0.697 0.687 0.678 0.668 0.658
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.648 0.638 0.628 0.618 0.608 0.599 0.589 0.579 0.569 0.559 0.548 0.538 0.528
##
## [[12]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.974 0.949 0.930 0.910 0.886 0.860 0.836 0.814 0.795 0.778 0.764 0.750
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.737 0.725 0.714 0.702 0.691 0.679 0.668 0.657 0.646 0.636 0.627 0.618 0.609
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.601 0.593 0.586 0.578 0.570 0.561 0.553 0.544 0.534 0.524 0.514 0.503 0.493
##
## [[13]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.972 0.941 0.907 0.865 0.819 0.774 0.735 0.702 0.675 0.654 0.638 0.625
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.615 0.608 0.601 0.595 0.590 0.585 0.580 0.576 0.572 0.569 0.567 0.565 0.563
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.561 0.559 0.556 0.551 0.545 0.536 0.527 0.515 0.503 0.490 0.477 0.464 0.451
##
## [[14]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.967 0.923 0.877 0.827 0.776 0.731 0.695 0.669 0.651 0.640 0.634 0.633
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.633 0.634 0.634 0.632 0.629 0.625 0.620 0.615 0.609 0.603 0.598 0.592 0.585
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.577 0.569 0.560 0.549 0.538 0.525 0.513 0.500 0.487 0.475 0.463 0.452 0.441
##
## [[15]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.974 0.951 0.932 0.914 0.895 0.875 0.853 0.831 0.809 0.789 0.770 0.752
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.735 0.719 0.704 0.689 0.675 0.661 0.649 0.636 0.624 0.612 0.601 0.589 0.577
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.566 0.555 0.544 0.533 0.523 0.513 0.504 0.495 0.486 0.477 0.469 0.461 0.453
##
## [[16]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.977 0.953 0.931 0.911 0.890 0.868 0.847 0.826 0.806 0.788 0.771 0.756
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.741 0.728 0.715 0.703 0.692 0.682 0.672 0.662 0.652 0.642 0.631 0.622 0.612
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.602 0.592 0.583 0.574 0.566 0.557 0.549 0.540 0.532 0.524 0.516 0.508 0.501
##
## [[17]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.970 0.933 0.894 0.846 0.793 0.739 0.688 0.642 0.603 0.569 0.543 0.523
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.510 0.503 0.501 0.503 0.509 0.517 0.525 0.533 0.539 0.542 0.543 0.540 0.535
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.527 0.517 0.505 0.493 0.480 0.467 0.455 0.443 0.432 0.421 0.409 0.398 0.387
##
## [[18]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.948 0.842 0.712 0.580 0.464 0.370 0.298 0.246 0.212 0.192 0.183 0.185
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.196 0.215 0.235 0.254 0.267 0.274 0.277 0.279 0.281 0.286 0.292 0.300 0.308
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.314 0.317 0.317 0.313 0.307 0.300 0.291 0.281 0.269 0.255 0.241 0.228 0.218
##
## [[19]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.985 0.965 0.941 0.916 0.889 0.860 0.832 0.805 0.781 0.761 0.745 0.731
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.721 0.713 0.706 0.701 0.696 0.692 0.687 0.683 0.678 0.672 0.666 0.660 0.653
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.646 0.639 0.631 0.624 0.616 0.608 0.600 0.592 0.584 0.575 0.566 0.557 0.547
##
## [[20]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.986 0.967 0.946 0.923 0.900 0.876 0.852 0.830 0.810 0.792 0.777 0.764
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.753 0.743 0.734 0.726 0.718 0.711 0.704 0.697 0.690 0.683 0.675 0.667 0.658
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.649 0.639 0.630 0.620 0.610 0.601 0.592 0.583 0.574 0.566 0.557 0.549 0.540
##
## [[21]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.986 0.969 0.952 0.935 0.916 0.897 0.878 0.858 0.839 0.821 0.803 0.787
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.771 0.757 0.744 0.731 0.719 0.708 0.697 0.687 0.677 0.668 0.658 0.649 0.640
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.631 0.622 0.613 0.604 0.595 0.586 0.577 0.568 0.558 0.549 0.540 0.531 0.521
##
## [[22]]
##
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1.000 0.986 0.971 0.955 0.940 0.924 0.908 0.892 0.875 0.860 0.846 0.833 0.820
## 13 14 15 16 17 18 19 20 21 22 23 24 25
## 0.809 0.798 0.787 0.776 0.765 0.755 0.745 0.734 0.724 0.714 0.704 0.695 0.685
## 26 27 28 29 30 31 32 33 34 35 36 37 38
## 0.676 0.667 0.658 0.650 0.641 0.632 0.624 0.615 0.606 0.596 0.587 0.577 0.567